1 Чтение и предобработка данных

Подгружаем датасет

insurance <- read_csv("./insurance_cost.csv")
## Rows: 1338 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): sex, smoker, region
## dbl (4): age, bmi, children, charges
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Ради большего интереса переименуем его.

endurance <- insurance

What endurance looks like

Смотрим на переменные

str(endurance)
## spc_tbl_ [1,338 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ age     : num [1:1338] 19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : chr [1:1338] "female" "male" "male" "male" ...
##  $ bmi     : num [1:1338] 27.9 33.8 33 22.7 28.9 ...
##  $ children: num [1:1338] 0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : chr [1:1338] "yes" "no" "no" "no" ...
##  $ region  : chr [1:1338] "southwest" "southeast" "southeast" "northwest" ...
##  $ charges : num [1:1338] 16885 1726 4449 21984 3867 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   age = col_double(),
##   ..   sex = col_character(),
##   ..   bmi = col_double(),
##   ..   children = col_double(),
##   ..   smoker = col_character(),
##   ..   region = col_character(),
##   ..   charges = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>
summary(endurance)
##       age            sex                 bmi           children    
##  Min.   :18.00   Length:1338        Min.   :15.96   Min.   :0.000  
##  1st Qu.:27.00   Class :character   1st Qu.:26.30   1st Qu.:0.000  
##  Median :39.00   Mode  :character   Median :30.40   Median :1.000  
##  Mean   :39.21                      Mean   :30.66   Mean   :1.095  
##  3rd Qu.:51.00                      3rd Qu.:34.69   3rd Qu.:2.000  
##  Max.   :64.00                      Max.   :53.13   Max.   :5.000  
##     smoker             region             charges     
##  Length:1338        Length:1338        Min.   : 1122  
##  Class :character   Class :character   1st Qu.: 4740  
##  Mode  :character   Mode  :character   Median : 9382  
##                                        Mean   :13270  
##                                        3rd Qu.:16640  
##                                        Max.   :63770

Все выглядит пристойно.

Оценим пропущенные значения.

endurance %>% 
  is.na() %>% 
  sum()
## [1] 0

Их нет, и это прекрасно.

EDA

2 Строим гистограммы

# Функция для массового создания гистограм в чанке setup
endurance_hsits <- lapply(endurance %>% select(where(is.numeric)) %>% colnames, function(x) histCreator(endurance, x))
endurance_hsits

  1. Распределение возраста можно было бы назвать равномерным, если бы не пики в районе 45-ти и 20-ти лет.
  2. Распределение ИМТ имееет колоколообразную форму с пиков районе 30-ти. Хвост в стороне высоких значений несколько больше.
  3. Распределение числа детей дискретное (логично!), убывающее на всем протяжении (были даже мысли перевести эту переменную в факторный вид, но я от них воздержался).
  4. Распределение суммы страховых выплат асимметричное, резко скошенное влево (пик в области низких значений), с длинным хвостом, уходящим в сторону больших значений.

3 Создаем график плотности

В данных не было указано, в какой валюте осуществлялись страховые выплаты. Выберем в качестве валюты наиболее вероятный вариант – монгольские тугрики.

charges_density <- ggplot(endurance)+
  geom_density(aes(x=charges),
               fill="white",
               color="green",
               alpha=0.9)+
  theme_dark()+   # Помню, что обычно рекомендуют theme_bw(), но раз в задании рекомендация поиграть с темами -- пусть будет соответствовать названию датасета
  geom_vline(aes(xintercept=mean(charges)),
             color="red",
             size=1.5,
             linetype="dashed")+
  geom_vline(aes(xintercept=median(charges)), 
             color="blue",
             size=1.5,
             linetype="dashed")+
 annotate(geom = "text", 
          x = 40000, 
          y = 0.00005, 
          label = glue("Среднее для страховых выплат составляет\n{round(mean(endurance$charges),0)} монгольских тугриков"),
          color = "red",
          size=5)+
   annotate(geom = "text", 
          x = 40000, 
          y = 0.00003, 
          label = glue("Медиана для страховых выплат составляет\n{round(median(endurance$charges),0)} монгольских тугриков"),
          color = "blue",
          size=5)+
  labs(x = "Страховые выплаты",
         y = "Плотность вероятности",
         title = "Рис. 2 График плотности вероятности\nстраховых выплат")+
    theme(axis.text.y = element_text(size=16),
        axis.title.x = element_text(size=19),
        axis.title.y = element_text(size=19),
        plot.title = element_text(size=20, hjust=0.5),
        legend.title = element_text(size=21),
        legend.text = element_text(size=18))
charges_density

4 Боксплоты

# Функция для массового создания боксплотов в чанке setup
endurance_boxps <- lapply(endurance %>% select(where(is.character)) %>% colnames, function(x) boxpCreator(endurance, x, "charges"))
endurance_boxps

5 Боксплоты и график плотности

boxplots <- ggarrange(endurance_boxps[[1]], 
          endurance_boxps[[2]],
          endurance_boxps[[3]],
          ncol=3, nrow=1)
combo_plot <- ggarrange(charges_density, boxplots,
          ncol=1, nrow=2)
annotate_figure(combo_plot, top = text_grob("Рис.4 Вероятность страховых выплат различного размера\nи распределение страховых выплат в зависимости от пола, курения и региона проживания\n",
               color = "black", face = "italic", size = 25))

6 Фасет по регионам

Базовый вариант

charges_density+
  facet_grid("region")+
  labs(title = "Рис. 5 График плотности страховых выплат по регионам")

7 Скаттерплот 1 (Возраст)

age_charges_scatter <- ggplot(endurance)+
  geom_point(aes(age, charges), colour="violet")+
  theme_bw()+
    labs(x = "Возраст",
         y = "Сумма страховки",
         title = "Рис. 6 Диаграмма рассеяния\nвозраст-страховые выплаты")+
    theme(axis.text.x = element_text(size=14),   # Как заказывали
        axis.text.y = element_text(size=10),
        axis.title.x = element_text(size=10), 
        axis.title.y = element_text(size=10),
        plot.title = element_text(size=20, hjust=0.5),
        legend.title = element_text(size=21),
        legend.text = element_text(size=18))
age_charges_scatter

8 Линия тренда

age_charges_scatter+
  geom_smooth(aes(x=age, y=charges), 
              color = "green",
              linewidth = 1.5,
              method = lm,
              se = FALSE)

9 Больше линий тренда!

age_charges_scatter+
  geom_point(aes(age, charges, colour=smoker))+
  geom_smooth(aes(age, charges, linewidth=smoker, group=smoker), 
              method=lm, 
              se=FALSE)+
  labs(colour="Курение",
       size="Курение")

10 Скаттерплот 2 (ИМТ)

ggplot(endurance, aes(bmi, charges, colour=smoker))+
  geom_point()+
  geom_smooth(aes(group=smoker), 
              method=lm, 
              se=FALSE)+
  theme_bw()+
    labs(x = "bmi",
         y = "Сумма страховки",
         title = "Рис. 7 Диаграмма рассеяния\nИМТ-страховые выплаты",
         colour="Курение")+
    theme(axis.text.x = element_text(size=18),
        axis.text.y = element_text(size=18),
        axis.title.x = element_text(size=20), 
        axis.title.y = element_text(size=20),
        plot.title = element_text(size=20, hjust=0.5),
        legend.title = element_text(size=21),
        legend.text = element_text(size=18))

Я бы сказал, что график по форме напоминает два шашлыка на шампуре под углом 35 градусов друг к другу. Из литературных данных известно, что подобная shashlyk-shaped форма графика наглядно показывает, что курение, вероятно, намного сильнее взаимосвязано с суммой страховых расходов, чем индекс массы тела.

Вопросы к данным

11 Вопрос №1.

Вопрос не про страховку: отличается ли ИМТ у имеющих одного ребенка некурящих мужчин и женщин.

Визуализируем данные

question1 <- endurance %>% 
  filter(children == 1 & smoker == "no") %>% 
  ggplot()+
    geom_boxplot(aes(x=sex, y=bmi, fill=sex))+
    theme_bw()+
    labs(x = "Пол",
         y = "ИМТ",
         title = "Рис. 7 Распределение ИМТ по полу у некурящих\nлиц, имеющих одного ребенка",
         fill = "Пол")+
    theme(axis.text.x = element_text(size=18),
        axis.text.y = element_text(size=18),
        axis.title.x = element_text(size=20), 
        axis.title.y = element_text(size=20),
        plot.title = element_text(size=20, hjust=0.5),
        legend.title = element_text(size=21),
        legend.text = element_text(size=18))
question1

Визуально кажется, что нет. Проверим.

endurance %>% 
  filter(children == 1 & smoker == "no") %>% 
  group_by(sex) %>% 
  do(broom::tidy(shapiro.test(.$bmi)))
## # A tibble: 2 × 4
## # Groups:   sex [2]
##   sex    statistic p.value method                     
##   <chr>      <dbl>   <dbl> <chr>                      
## 1 female     0.994  0.837  Shapiro-Wilk normality test
## 2 male       0.981  0.0661 Shapiro-Wilk normality test

Не можем отвергнуть нулевую гипотезу о том, что распределения ИМТ в целевых группах нормальны. Используем тест Стьюдента для проверки гипотезы о равенстве средних.

stat.test <- endurance %>% 
  filter(children == 1 & smoker == "no") %>% 
  t.test(bmi ~ sex, data = .) #%>% 
  #add_xy_position(x = "sex")
stat.test
## 
##  Welch Two Sample t-test
## 
## data:  bmi by sex
## t = -0.95426, df = 256.89, p-value = 0.3408
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -2.2035331  0.7650241
## sample estimates:
## mean in group female   mean in group male 
##             30.20936             30.92862

Нулевую гипотезу отвергнуть не можем. Средние по нашим данным не различаются.

Обоснование типа графика

Имеем одну категориальную и одну непрерывную числовую переменную. Оптимальным способом отображения данных такого типа являются графики типов боксплот, скрипкаграфик вайолинплот и джиттерплот

12 Вопрос №2.

Как связаны между собой возраст, размер страховых выплат и ИМТ у курящих и некурящих мужчин?

Визуализируем данные

plot_ly(
  endurance %>% filter(sex == "male"),
  x = ~ age,
  y = ~ bmi,
  z = ~ charges, 
  size = 1,
  color = ~ smoker,
  colors = c("#FF0033", "#00CC00"),
  alpha = 0.5,
  width = 1.5,
  height = 1.5
) %>% 
  layout(
      title = "Распределение возраста, ИМТи размера выплаченной страховки\nу курящих и некурящих мужчин",
      xaxis = list(title = "Возраст",
                   zeroline = FALSE),
      yaxis = list(title = "Индекс массы тела",
                   zeroline = FALSE)
  )

Вырисовываются три паралельные группы, в которых с увеличением возраста растет и сумма выплаченной страховки (как у курящих, так и у некурящих, причем у последних суммы были заметно ниже). В случае с индексом массы тела его увеличение приводило к увеличению страховки только у некурящих пациентов.

Обоснование типа графика

Имеем три типа переменных. Посмотреть на их взаимное распределение удобно на диаграмме рассеяния (скаттерплоте). Ну а визуализация в плотли выбрана потому, что это забавно выглядит (этот пункт был выполнен до публикации второй домашки, а после передлывать как-то не пришлось по душе).

13 Вопрос №3.

Как предсказать размер страховки, исходя из имеющихся данных?

Попробуем поиграть с линейной регрессией.

# Стандартизуем количественные переменные
endurance_scaled <- endurance %>% 
  mutate(across(where(is.numeric), scale))
# Построим модель
charge_lm <- lm(charges ~ bmi + age + children + sex + smoker + region, data = endurance_scaled)
summary(charge_lm)
## 
## Call:
## lm(formula = charges ~ bmi + age + children + sex + smoker + 
##     region, data = endurance_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9335 -0.2352 -0.0811  0.1151  2.4767 
## 
## Coefficients:
##                 Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)     -0.34822    0.03189 -10.920 < 0.0000000000000002 ***
## bmi              0.17081    0.01440  11.860 < 0.0000000000000002 ***
## age              0.29800    0.01380  21.587 < 0.0000000000000002 ***
## children         0.04733    0.01372   3.451             0.000577 ***
## sexmale         -0.01084    0.02749  -0.394             0.693348    
## smokeryes        1.96932    0.03412  57.723 < 0.0000000000000002 ***
## regionnorthwest -0.02915    0.03933  -0.741             0.458769    
## regionsoutheast -0.08547    0.03953  -2.162             0.030782 *  
## regionsouthwest -0.07928    0.03947  -2.009             0.044765 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5006 on 1329 degrees of freedom
## Multiple R-squared:  0.7509, Adjusted R-squared:  0.7494 
## F-statistic: 500.8 on 8 and 1329 DF,  p-value: < 0.00000000000000022

Выглядит неплохо уже на текущем этапе – поправленный R-квадрат = 0.75

Поищем мультиколлинеарность в данных (вот здесь мы и применим визуализацию)

Визуализируем данные

endurance_corr <- corr.test(endurance %>% select(where(is.numeric)& !charges), method="spearman")
corr_plot <- ggcorrplot(endurance_corr$r, 
                        type = "lower", 
                        outline.col = "white", 
                        lab=TRUE, 
                        method = "circle", 
                        p.mat = endurance_corr$p,
                        title = "Скоррелированность количественных данных\nв датасете по страховке",
                        legend.title = "Коэффициент\nкорреляции")
corr_plot+
  labs()

На первый взгляд коллинеарности нет. Проверим формальными тестами.

vif(charge_lm)
##              GVIF Df GVIF^(1/(2*Df))
## bmi      1.106630  1        1.051965
## age      1.016822  1        1.008376
## children 1.004011  1        1.002003
## sex      1.008900  1        1.004440
## smoker   1.012074  1        1.006019
## region   1.098893  3        1.015841

Все действительно в порядке.

Сделаем шаг backward selection:

drop1(charge_lm, test = "F")
## Single term deletions
## 
## Model:
## charges ~ bmi + age + children + sex + smoker + region
##          Df Sum of Sq     RSS      AIC   F value                Pr(>F)    
## <none>                 333.03 -1842.76                                    
## bmi       1     35.25  368.28 -1710.15  140.6627 < 0.00000000000000022 ***
## age       1    116.77  449.80 -1442.60  465.9837 < 0.00000000000000022 ***
## children  1      2.98  336.01 -1832.82   11.9063              0.000577 ***
## sex       1      0.04  333.07 -1844.60    0.1556              0.693348    
## smoker    1    834.95 1167.98  -165.84 3331.9680 < 0.00000000000000022 ***
## region    3      1.59  334.62 -1842.38    2.1173              0.096221 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Удалим пол, как наименее значимый предиктор.

charge_lm2 <- update(charge_lm, .~. - sex)
drop1(charge_lm2, test = "F")
## Single term deletions
## 
## Model:
## charges ~ bmi + age + children + smoker + region
##          Df Sum of Sq     RSS      AIC   F value                Pr(>F)    
## <none>                 333.07 -1844.60                                    
## bmi       1     35.22  368.28 -1712.12  140.6226 < 0.00000000000000022 ***
## age       1    116.95  450.02 -1443.95  466.9969 < 0.00000000000000022 ***
## children  1      2.97  336.04 -1834.71   11.8706              0.000588 ***
## smoker    1    838.82 1171.89  -163.37 3349.5461 < 0.00000000000000022 ***
## region    3      1.59  334.66 -1844.23    2.1166              0.096315 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Еще шаг. Теперь удаляем регион.

charge_lm3 <- update(charge_lm2, .~. - region)
drop1(charge_lm3, test = "F")
## Single term deletions
## 
## Model:
## charges ~ bmi + age + children + smoker
##          Df Sum of Sq     RSS      AIC  F value                Pr(>F)    
## <none>                 334.66 -1844.23                                   
## bmi       1     34.70  369.36 -1714.24  138.203 < 0.00000000000000022 ***
## age       1    117.94  452.60 -1442.28  469.789 < 0.00000000000000022 ***
## children  1      2.96  337.62 -1834.43   11.809             0.0006077 ***
## smoker    1    841.77 1176.43  -164.19 3352.911 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(charge_lm3)
## 
## Call:
## lm(formula = charges ~ bmi + age + children + smoker, data = endurance_scaled)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.98249 -0.24119 -0.08147  0.11497  2.43679 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) -0.40266    0.01536 -26.211 < 0.0000000000000002 ***
## bmi          0.16207    0.01379  11.756 < 0.0000000000000002 ***
## age          0.29916    0.01380  21.675 < 0.0000000000000002 ***
## children     0.04713    0.01372   3.436             0.000608 ***
## smokeryes    1.96626    0.03396  57.904 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5011 on 1333 degrees of freedom
## Multiple R-squared:  0.7497, Adjusted R-squared:  0.7489 
## F-statistic: 998.1 on 4 and 1333 DF,  p-value: < 0.00000000000000022

Наш R-квадрат почти не изменился, зато теперь нам достаточно четырех переменных вместо шести (что весьма важно – мы исключили достаточно трудно интерпретируемую в рамках линейной регрессию переменную region).

Обоснование типа графика

Коррплот сделан как раз для того, чтобы искать скоррелированные переменные. В нашем случае его наглядность может несколько страдать из-за малого количества числовых переменных, одако, когда у нас 10-15 оных, коррплоты могут очень облегчить поиск коллинеарных переменных.

14 Воспроизведение графика

ВАЖНО

У целевого графика есть расхождения с подписями. Если делать его с разбиением на группы в соответствии с легендой, то точек на скаттерплоте меньше, чем на целевом. Если же сбросить нижнюю границу возраста в первой группе (что приведет к включению лиц в возрасте 18-20 лет в группу 21-34), график принимает необходимый вид. Поскольку целью задания является именно построение идентичного графика, я занулил нижнюю границу в первой возрастной группе.

График

endurance %>% 
  mutate(age_group=as.factor(case_when(
    age < 35 ~ "age: 21-34",
    age >= 35 & age < 50 ~ "age: 35-49",
    age >= 50 ~ "age: 50+"
  ))) %>% 
  na.omit() %>% 
  ggplot(aes(x=bmi, y=log(charges)))+
    geom_point(color = "#330066", alpha = 0.5, size=2)+
    facet_grid(cols=vars(age_group))+
    geom_smooth(aes(colour=age_group), method=lm)+
  theme_minimal()+
  labs(title="Отношение индекса массы тела к логарифму трат по возрастным группам")+
  theme(
    axis.text.x = element_text(size=16),
    axis.text.y = element_text(size=16),
    axis.title.x = element_text(size=19), 
    axis.title.y = element_text(size=19),
    plot.title = element_text(size=23, hjust=5),
    legend.title = element_text(size=19),
    legend.text = element_text(size=16),
    legend.position = "bottom",
    legend.key.size = unit(1.5, "line"),
    strip.text.x = element_text(size=16)
  )

Смею тешить себя надеждой, что